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Abstract 

We have performed direct numerical simulation of turbulent open channel flow over a 
smooth horizontal wall in the presence of finite-size, heavy particles. The spherical particles 
have a diameter of approximately 7 wall units, a density of 1.7 times the fluid density 
and a solid volume fraction of 5 ■ 10 -4 . The value of the Galileo number is set to 16.5 
, while the Shields parameter measures approximately 0.2. Under these conditions, the 
particles are predominantly located in the vicinity of the bottom wall, where they exhibit 
strong preferential concentration which we quantify by means of Voronoi analysis and by 
computing the particle-conditioned concentration field. As observed in previous studies 
with similar parameter values, the mean streamwise particle velocity is smaller than that 
of the fluid. We propose a new definition of the fluid velocity "seen" by finite-size particles 
based on an average over a spherical surface segment, from which we deduce in the present 
case that the particles are instantaneously lagging the fluid only by a small amount. The 
particle-conditioned fluid velocity field shows that the particles preferentially reside in the 
low-speed streaks, leading to the observed apparent lag. Finally, a vortex eduction study 
reveals that spanwise particle motion is significantly correlated with the presence of vortices 
with the corresponding sense of rotation which are located in the immediate vicinity of the 
near- wall particles. 

1 Introduction 

The transport of solid particles in wall bounded turbulent flows is a common occurrence in various 
natural and man-made systems. This transport phenomenon has significant implications in many 
environmental and industrial processes. For instance one of the key environmental variables in 
river geomorphology is sediment transport, which involves the erosion, movement and deposition 
of sediment particles by the fluid flow. An improved understanding of the coupled interaction 
between solid particles and turbulence is highly desirable as it would pave the way for improve- 
ments of engineering-type formulas. However the complex structure of wall-bounded shear flows 
and the dependence on multiple governing parameters has made this task very challenging to 
the present date. 

Experimental studies of particle-turbulence interaction in horizontal channel flow show that 
near-wall coherent structures play an important role in the dynamics of particle motion (see 
e.g. Sumer and Oguz, 1978; Yung, Merry, and Bott, 1989; Rashidi, Hetsroni, and Banerjee, 
1990; Kaftori, Hetsroni, and Banerjee, 1995a, b; Nino and Garcia, 1996; Kiger and Pan, 2002; 
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Figure 1: Geometrical configuration of the present simulation. Gravity is directed in the negative 
y- direction. The computational domain is periodic along the x- and z-directions. 

Righetti and Romano, 2004). Rashidi et al. (1990) investigated the coupled interaction between 
phases in a horizontal channel flow laden with heavy spherical particles. They observed that 
most of the particles, once they settle to the bottom wall, accumulate in the low-speed streaks. 
The particles, depending on their size, their density and flow Reynolds number, are observed 
to get lifted up and entrained into the outer flow presumably by the action of the coherent 
structures. Particles resettling to the bottom are observed to migrate into the low-speed fluid 
regions, supposedly by the action of the eddy structures. Hetsroni and Rozcnblit (1994) found 
similar particle segregation at the wall when investigating particles with a diameter of 10 wall 
units (i.e. equal to 10 times the viscous length scale v/u T where v is the kinematic viscosity 
of the fluid and u T is the friction velocity). However, larger particles having a diameter of 
more than 30 wall units did not accumulate in the low-speed fluid regions; rather they formed a 
random distribution at the bottom wall of their horizontal flume. In similar experiments, Kaftori 
et al. (1995a) and Nino and Garcia (1996) investigated the effect of the coherent structures on 
particle motion, distribution, entrainment as well as deposition in the wall region. They show 
that heavy particles directly interact with the flow structures. At the bottom wall, particles were 
observed to be non-uniformly distributed and to form streamwise aligned streaks. The shape, 
length and persistence with time of these streaks was observed to vary with the particle size 
and/or the flow rate. Formation of particle streaks was also observed by Yung et al. (1989) 
even for particles which are completely submerged within the viscous sublayer (diameter smaller 
than 1.3 wall units). The experimental studies of Kaftori et al. (1995b); Kigcr and Pan (2002); 
Righetti and Romano (2004) were focused on the characteristics of suspended heavy particles and 
their turbulence modulation effects in open as well as closed channel flows. These experiments 
showed that generally, the mean streamwise velocity of the particles is smaller than that of the 
fluid except for those particles located in a layer very close to the wall where contrarily particle 
velocity is reported to be on average higher than that of the fluid. 

The majority of available numerical studies on particle-turbulence interaction are based on 
the point-particle approach (see e.g. Rouson and Eaton, 2001; Marchioli and Soldati, 2002; 
Narayanan, Lakehal, Botto, and Soldati, 2003; Picciotto, Marchioli, and Soldati, 2005). Studies 
based on this approach show that there is a strong correlation between the near-wall coherent 
structures and particle deposition and entrainment processes as well as the resulting preferen- 
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Table 1: Physical parameters of the simulation. Re^ is the Reynolds number based on the 
channel depth h and bulk velocity Mb- Re T is the Reynolds number based on h and the friction 
velocity u T . St + and St^ are the Stokes numbers based on near wall fluid time scale and bulk 
time scale respectively. D is the diameter of the spherical particles. $ s is the global solid volume 
fraction and |g| is the magnitude of acceleration due to gravity. Ga and 6 are the Galileo number 
and Shields number, respectively. 

tial particle segregation. The degree of the correlation is found to be dependent on the Stokes 
number, i.e. the ratio of the particle response time to some representative time scale of the 
flow (Soldati and Marchioli, 2009). However, the point-particle approach is limited to particle 
sizes which are much smaller than the smallest flow scales. Therefore, in order to numerically 
investigate finite-size effects, one has to resort to simulations which fully resolve the particles. 

Recently more faithful simulations in which the interface between the phases is resolved 
to accurately compute the flow field at the particle scale are starting to emerge both for wall 
bounded flows (Pan and Banerjee, 1997; Uhlmann, 2008; Shao, Wu, and Yu, 2012; Garci'a- 
Villalba, Kidancmariam, and Uhlmann, 2012) and unbounded flows (Kajishima and Takiguchi, 
2002; Ten Cate, Derksen, Portcla, and Van Den Akker, 2004; Lucci, Ferrante, and Elghobashi, 
2010, 2011). In one of the first simulations of this kind, Pan and Banerjee (1997) simulated 
horizontal open channel flow seeded with a small number of fixed and moving particles of diameter 
8.5 and 17 wall units at an average solid volume fraction of O(10 -4 ). They observed that the 
presence of particles strongly affected the turbulent structures. Uhlmann (2008) and Garcia- 
Villalba ct al. (2012) simulated vertically oriented channel flow seeded with freely moving heavy 
particles with a diameter of approximately 11 wall units at a solid volume fraction of 0.4%. These 
simulations revealed that the particles strongly modified the flow and lead to the formation of 
very large streamwise-elongated flow structures. Recently, Shao et al. (2012) simulated particle- 
laden turbulent flow in a horizontal channel. They considered heavy particles (with density 1.5 
times that of the fluid) which are relatively large in size (with diameter 10% and 20% of the 
half channel height), at relatively large solid volume fraction values (up to 7%). They reported 
that, in cases where settling effects are not considered (without gravity), the particle distribution 
was homogeneous and the mean fluid and particle velocities were equal. However, when settling 
effects are considered, they reported that most particles settle to the bottom wall where they 
accumulate in the low-speed fluid regions. 

The general picture that has emerged from previous studies of horizontal turbulent channel 
flow with heavy particles (in the regime where streaky particle patterns are observed) can be 
summed up briefly as follows. Sweeps of high-speed fluid as well as gravitational settling bring 
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Table 2: Numerical parameters of the simulation. Lj and Ni are the computational domain 
length and number of grid points employed in the ith coordinate direction, respectively. N p is 
the number of particles considered. T b s is the observation time of the simulations excluding the 
transient period. 
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Figure 2: Time evolution of (a) box-averaged turbulent kinetic energy and (b) friction- velocity 
based Reynolds number. Note that the time coordinate has been shifted to the beginning of the 
observation interval, excluding initial transients. 

suspended heavy particles towards the wall. Through the action of the quasi-streamwise vor- 
tices, particles near the wall are forced to move in the spanwise direction into the low-speed 
regions, resulting in particle accumulation in the form of persistent streaky particle patterns. 
The counter-rotating quasi-streamwise vortices flanking the low-speed streaks generate ejections 
of low-momentum fluid away from the wall. These ejection events are the predominant sources of 
upward particle motion from the wall. Concerning the difference between the mean fluid velocity 
and the mean particle velocity (i.e. particles apparently lagging the fluid), it is believed to be a 
statistical consequence of particles preferentially residing in the low-speed regions, rather than a 
manifestation of particles instantaneously lagging the surrounding fluid in a significant manner. 
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Despite the great progress achieved through past studies, a number of important questions 
still remain unanswered. The present work is an attempt to fill these gaps. First, a detailed 
characterization of the spatial structure of the dispersed phase in horizontal channel flow with 
finite-size particles is still lacking. For this purpose an analysis based upon Voronoi tesselation 
relative to the instantaneous particle distribution is an adequate instrument, as has been demon- 
strated by Monchaux, Bourgoin, and Cartellier (2010) and Monchaux, Bourgoin, and Cartellier 
(20f2). Additionally, the particle-pair distribution function yields complementary information 
on the large-scale structure of the particle 'field'. Second, in order to rigorously confirm the 
explanation of the apparent velocity lag of the particles in terms of preferential sampling of the 
fluid velocity field, the possibility of a systematic instantaneous inter-phase velocity lag needs to 
be excluded. This in turn requires the definition of a characteristic fluid velocity in the vicinity 
of the particles (often referred to as the velocity 'seen' by the particles). However, for finite-size 
particles a unique definition does not exist, and for different reasons previous attempts (Bagchi 
and Balachandar, 2003; Merle, Legendre, and Magnaudct, 2005; Zeng, Balachandar, Fischer, and 
Najjar, 2008; Lucci et al., 2010) do not yield the desired results in the present flow configuration. 
Here we propose a new definition of the fluid velocity seen by finite-size particles. A third point 
of interest concerns the fluid velocity field conditioned upon the presence of particles and, in 
particular, the relation between the spanwise motion of near- wall particles (driving particles into 
low-speed streaks) and the presence of coherent vortices. 

In order to investigate these three aspects of the problem, highly-resolved flow and particle 
data (in time and space) is required. In the present work we have generated such data by means 
of interface-resolved direct numerical simulation of particle-laden, horizontal, open channel flow 
over a smooth wall. All relevant scales of the flow problem are resolved by means of a finite- 
difference/immersed-boundary technique. The solid volume fraction was set to a relatively low 
value in order to avoid a dominance of inter-particle collisions. Contact between pairs of particles 
and contact of particles with the solid wall is considered as frictionless in the simulations. 

2 Computational methodology and setup 
2.1 Numerical method 

The numerical method employed in the present study is a formulation of the immersed boundary 
method for the simulation of particulate flows developed by Uhlmann (2005a). The basic idea of 
the immersed boundary method is to solve the modified incompressible Navier-Stokes equations 
throughout the entire domain Q comprising the fluid domain Qf and the domain occupied by the 
suspended particles tt s while adding a force term which serves to impose the no-slip condition 
at the fluid-solid interface. The immersed boundary technique is realized in the framework of a 
standard fractional step method for the incompressible Navier-Stokes equations. The temporal 
discretization is semi-implicit, based on the Crank-Nicolson scheme for the viscous terms and 
a low-storage three-step Runge-Kutta procedure for the non-linear part (Verzicco and Orlandi, 
1996). The spatial operators are evaluated by central finite-differences on a staggered grid. 
The temporal and spatial accuracy of this scheme are of second order. In the computation of 
the forcing term, the necessary interpolation of variable values from Eulerian grid positions to 
particle-related Lagrangian positions (and the inverse operation of spreading the computed force 
terms back to the Eulerian grid) are performed by means of the regularized delta function ap- 
proach of Pcskin (2002). This procedure yields a smooth temporal variation of the hydrodynamic 
forces acting on individual particles while they are in arbitrary motion with respect to the fixed 
grid. The employed Cartesian grid is uniform and isotropic in order to ensure the conservation of 
important quantities such as the total force and torque, during the interpolation and spreading 
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Figure 3: Instantaneous three-dimensional snapshot of the flow field and the particle positions 
(observed mostly residing at the bottom wall). Strong vortical regions are shown by plotting a 
negative- valued iso-surface of the A2 field (cf. definition in § 3.6). The iso-surfaces are colored 
(red/blue) based on the sign of the streamwise vorticity (positive/negative). The greyscales 
indicate values of the streamwise velocity fluctuation in a wall-parallel plane at y + = 5 (projected 
upon the plane y = for better visibility) and in wall perpendicular planes at x/h = 12 and 
z/h = 0, ranging from black (minimum negative values) to white (maximum positive values). 
Note that although the particle surfaces are colored in green, particles appear mostly in either 
red or blue color, due to vortical structures located in their immediate vicinity. A close-up of 
the same graph looking against the spanwise direction is shown in (b) (i.e. the flow is from left 
to right). 

procedure. The particle motion is determined by the Runge-Kutta-discretized Newton equations 
for linear and angular motion of rigid bodies, driven by buoyancy, hydrodynamic forces/torque 
and contact forces (in case of collisions) . 

During the course of a simulation, particles can closely approach each other. However, very 
thin inter-particle fluid films cannot be resolved by a typical grid, and therefore, the correct 
buildup of repulsive pressure is not captured, which, in turn, can lead to possible partial "overlap" 
of the particle positions in the numerical computation. In practice, particle-particle contacts are 
treated by a simple repulsive force mechanism (Glowinski, Pan, Hesla, and Joseph, 1999), which 
is active for particle pairs at a distance smaller than two grid spacings. The contact force only 
acts in the normal direction (along the line connecting the two particles' centers) and frictional 
contact forces are not considered. The analogous treatment is applied to particle-wall contact. 
Note that a posteriori evaluation of the particle trajectories revealed that the average temporal 
interval between two collision events measures approximately 4.1 bulk flow time units. 

The numerical method employs domain decomposition for parallelism and has been shown 
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Figure 4: Profile of the average solid volume fraction (<fi s ), normalized by the global value <& s . 
The width of the averaging bins are shown with gray-colored lines. The number of data samples 
located in bins above the chain-dotted line (at y + = 36) is smaller than 800. Note that the bin 
adjacent to the wall has zero number of samples as there is no particle center in the interval 
< y P < D/2. The insert shows the same quantity for an interval close to the wall, evaluated 
with finer averaging bins Ah = Ax. 

to run on grids of up to 8192 3 , using up to O(10 5 ) processor cores in scaling tests (Uhlmann, 
2010). The numerical method has been validated on a whole range of benchmark problems 
(Uhlmann, 2004, 2005a,b, 2006), and has been previously employed for the simulation of several 
flow configurations (Uhlmann, 2008; Chan-Braun, Garcia- Villalba, and Uhlmann, 2011; Garcfa- 
Villalba et al., 2012). 

2.2 Flow configuration and parameter values 

Horizontal open channel flow laden with finite-size, heavy, spherical particles is considered in- 
cluding the action of gravity. As shown in figure 1, a Cartesian coordinate system is adopted 
such that x, y, and z are the streamwise, wall-normal and spanwise directions respectively. Mean 
flow and gravity are directed in the positive x and the negative y directions respectively. The 
computational domain is periodic in the homogeneous directions (streamwise and spanwise). A 
free-slip condition is imposed at the top boundary while the no-slip condition is imposed at the 
bottom wall. The channel is driven by a horizontal mean pressure gradient imposing a constant 
flow rate. The Reynolds number of the flow based on the channel height h and the bulk velocity 
Mb = r J^(u[)dy is set to a value Re^, = u^h/v = 2870 generating a turbulent flow at a friction 
Reynolds number Re T = u T hjv — 184, where u T = yj (r w ) / p is the friction velocity, and (t w ) 
is the mean shear stress averaged in time and over the solid wall at y = 0. The value of the 
friction velocity matches the value obtained in single-phase flow at the same bulk Reynolds num- 
ber to within 1%. The spherical particles considered are relatively small in size when compared 
to the channel depth (particle diameter D/h = 10/256) but they are larger than the viscous 
length scale (D + = 7.21). Note that the " + " superscript indicates scaling in wall units. The 
particle-to-fluid density ratio is set to a value of p p / pf — 1.7. The Stokes number is defined as 
the ratio between the particle relaxation time scale t p = p p / pfD 2 / '(18^), assuming Stokes drag 
law, and a fluid time scale. The Stokes number based on the near wall fluid time scale has a 
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Figure 5: Wall-normal profile of the mean streamwise fluid velocity (solid line) and particle 
velocity (filled symbols) both in (a) outer units and (b) wall units. The dashed line represents 
the corresponding profile of a single phase flow at a bulk Reynolds number identical to the current 
simulation. 

value St + = T v v? T jv = 4.9 and the Stokes number based on bulk fluid time scale Tb = h/u^ has 
a value of S% = T p /Tb = 0.41. The global solid volume fraction was set to $ s = 5 x 1CP 4 . This 
means that there are N p = 518 particles in the computational domain fl of size 12h x h x 3/i. 
The non-dimensional gravity constant gh/u^ was set to a value 0.8 such that the Galileo number, 
defined as Ga 2 = (p P /pf — \)gD i jv 2 has a value Ga = 16.49 and the Shields number, defined 
as O^ 1 = (pp/pf — l)gD/u 2 has a value 9 = 0.19. The computational domain is discretized by 
a uniform isotropic grid with grid spacing Ax = Ay = Az = h/256. This grid spacing yields a 
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particle resolution of D/Ax = 10 and in terms of wall units Ax + ~ 0.7. Table 1 and table 2 
summarize the important physical and numerical parameter values adopted. 

The simulation was first initiated on a coarser grid (by a factor of two in each coordinate direc- 
tion), using a fully turbulent flow field with particles distributed throughout the computational 
domain as initial condition. A transient of roughly 40 Tb was observed, after which the particles 
had formed patterns very similar to those observed and described below. Then the final data 
of the coarser-grid simulation was interpolated upon the current grid, and the simulation was 
continued for another 252 Tb. The entire initial interval of 560 Tb was unintentionally computed 
with an incorrect value of the particles' moment of inertia (p p D 5 n5/A8 instead of p p D 5 'it /60). 
After correcting the value of the moment of inertia, the simulation was run for additional 260 Tb, 
the first 5 Tb of which have been discarded. Therefore, the temporal observation interval over 
which the present data was generated measures T D b s = 255 T>. 

Figure 2 shows the time evolution over the observation interval of the box-averaged turbulent 
kinetic energy, defined as (fc) xyz = §((itf,<itf,i)xyz - (wf,i) X yz(wf,i) xy z), where (-) xyz denotes an 
instantaneous average over the domain occupied by the fluid, as defined in (22). The time 
evolution of Re T based upon the instantaneous value of the friction velocity is also presented in 
the figure. Several cycles of large-scale fluctuations can be observed. The reader is referred to A 
for a statistical convergence check based on the balance of the mean streamwise momentum. 

Figure 3 gives an impression of the size of the domain, the coherent flow structures and the 
particle distribution. 

The simulation was typically performed on 48 to 256 processor cores and required roughly 
half a million core hours of CPU time. 

In addition, we have generated single-phase reference data at the same parameter point by 
means of a pseudo-spectral method (Kim, Moin, and Moser, 1987). This reference simulation has 
been run for over 4400 bulk time units, i.e. a considerably longer interval than that of the main 
simulation. The single-phase data will be used for the purpose of comparison in the following. 

2.3 Notation 

Before turning to the results, let us fix the basic notation followed throughout the present text. 
Velocity vectors and their components corresponding to the fluid and the particle phases are dis- 
tinguished by subscript "f" and "p" , respectively, as in Uf = (uf , Vf , Wf) T and u p = (u p , v p , w p ) T . 
Similarly, the particle position vector is denoted as x p = (x p ,y p , z p ) T and the vector of angular 
particle velocity as ui p = (oj px ,u! py ,u> PtZ ) T . 

Fluctuations of the fluid velocity field with respect to the average over wall-parallel planes 
and time are henceforth denoted by a single prime, i.e. Uf(x, t) = Uf(x, t) — (u[)(y). Likewise, 
the fluctuations of the particle velocity are defined as the difference between the instantaneous 
value and the average (over time and all particles located in predefined wall-normal intervals) at 
the corresponding location, viz. u p (x p (t),t) = u p (x p (i),i) — (u p )(y p (t)). The reader is referred 
to B for definitions of the various averaging operators used in the present study. 

The particle radius is henceforth denoted as R = D/2. 
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Figure 6: Wall-normal profile of the rms velocity fluctuations of the fluid phase (lines) and 
particles (filled symbols) in (a) outer units and (&) wall units. The dashed line indicates the 
single phase data. The black, blue and red colors represent the streamwise, spanwise and wall- 
normal components of the fluctuations, respectively. The insert shows the same quantities but 
the particle data is evaluated at finer bins of width Ah = Ax. 

3 Results 

3.1 Eulerian statistics 

3.1.1 Mean solid volume fraction 

The global solid volume fraction has such a low value that even if all particles are in contact 
with the solid wall, their projected area covers less than two percent of the wall-plane, i.e. 
NpirR 2 / (L X L Z ) = 0.0172. Figure 4 shows the wall-normal profile of the average solid volume 
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fraction (<p s ) (cf. averaging operator denned in equation 25 in B.3). Due to the gravitational 
settling effect, a strong wall- normal particle concentration gradient forms near the wall. The wall- 
normal location of the particle centers of the overwhelming majority of all particles (99.15%) is 
in the interval D/2 < y p < D (i.e. corresponding to the second averaging bin). At larger wall- 
distances, the number of particle samples becomes scarce for the given temporal observation 
interval. In order to ensure an adequate quality of the statistics, we have chosen to consider only 
those bins with y+ < 36 (up to the dashed line in figure 4) , for which the number of samples per 
bin ranges from a maximum of 10 7 down to 1386. Therefore, in the following presentation, all 
Eulerian particle statistics are only provided for wall distances y + < 36. 

3.1.2 Mean streamwise velocity 

The profiles of the mean streamwise velocity of fluid and solid phase are shown in figure 5 in 
outer and inner scaling. Additionally, the data for single-phase flow is included for comparison. 
In both scalings, the velocity profile of the fluid phase, (itf), almost coincides with the results of 
the single phase flow. The difference between the curves is smaller than 9.3 x 10~ 3 times the 
maximum velocity, and, therefore, within the range of statistical uncertainty. Thus, for the given 
parameters, the presence of particles has a negligible effect on the mean fluid velocity profile. 

The mean streamwise particle velocity, (u p ), is observed to be systematically smaller than 
that of the fluid phase (ttf). The difference in the mean velocities of the two phases, i.e. the 
apparent velocity lag, is denoted as 

"lag = (Uf) - (U p ) , (1) 

which is found to vary in the range of 2.1u T to 3.2u r , corresponding to 0.14ttb to 0.21ub (cf. 
figure 14& which will be discussed below). The general feature of a positive apparent velocity 
lag ttiag > has been reported in a number of experimental studies on horizontal wall-bounded 
shear flow (Rashidi ct al., 1990; Kaftori ct al., 1995b; Tanicre, Oesterle, and Monnier, 1997; Kiger 
and Pan, 2002; Righetti and Romano, 2004; Muste, Yu, Fujita, and Ettema, 2009; Noguchi and 
Nezu, 2009) as well as in the interface-resolved DNS of horizontal plane channel flow by Shao 
ct al. (2012). The aforementioned references cover a relatively broad range of parameter values in 
terms of Stokes number, Galileo number and solid volume fraction. Therefore, the phenomenon, 
which is commonly linked to particles residing preferentially in low-speed regions of the flow, 
appears to be a relatively 'robust' feature of particulate flow in horizontal, wall-bounded, shear 
turbulence. 

However, it should be mentioned that in most of the aforementioned experimental studies 
a change of sign of the apparent velocity lag is observed - contrary to the present simulations 
and contrary to those by Shao et al. (2012). Kaftori ct al. (1995b) measured particle velocities 
exceeding the average fluid velocity in a small interval (y < D) very near the wall. The authors, 
however, attribute this result to measurement errors (influence of the angular particle velocity) 
rather than to a physical effect. Kiger and Pan (2002) observed average particle velocities leading 
the average fluid velocity at wall-distances y + < 35; however, they did not further comment on 
the possible cause for this result. Righetti and Romano (2004) also detected a change of sign 
of the apparent slip velocity, i.e. (u p ) being larger than (itf) for y + < 20 (slightly depending on 
particle diameter), and vice versa further away from the wall. The authors explained the leading 
particle velocity near the wall as a consequence of sweep events (i.e. in the fourth quadrant of 
the Mp,Wp-plane) dominating the particle transport from the outer flow towards the wall. In the 
present case, on the contrary, the vast majority of the particles are residing in the direct vicinity 
of the wall (cf. concentration profile in figure 4) , and excursions of individual particles to larger 
wall-distances occur only very infrequently. Therefore, the importance of particles carrying high 



11 



(a) 




Figure 7: Wall- normal profile of the Reynolds shear stress of the fluid phase (lines) and correlation 
between streamwise and wall-normal velocity fluctuations of particles (filled symbols) in (a) outer 
units and (6) wall units. The dashed line indicates the single phase Reynolds stress profiles. The 
insert shows the same quantities but the particle data is evaluated at finer bins of width Ah = Ax. 

streamwise momentum while being swept towards the wall is of little importance to the mean 
streamwise velocity budget in the region near the wall. As a consequence, negative values of iti ag 
are not observed in the present case. 

In § 3.2 we will return to the apparent velocity lag in order to determine the contribution 
from the instantaneous lag with respect to a characteristic fluid velocity in the vicinity of the 
particles. In § 3.4 it will be shown that indeed the decisive contribution stems from particles 
preferentially sampling the low-speed regions of the flow. 
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Figure 8: Wall-normal profile of the number of collision events happening mainly in the vertical 
direction, i.e. the orientation of the line joining the centers of the two colliding particles is closer 
to a vertical line than it is to a horizontal line. The collision event is attributed to the particle 
which is at a larger wall-normal distance than its partner. Data points correspond to the same 
fine bins adopted in the insert of figure 7. 



3.1.3 Covariances of fluid and particle velocity fluctuations 

Figure 6 and figure 7 show the wall-normal profiles of the rms velocity fluctuation components 
of the fluid phase (u' f a u' f a ) 5 and of the particle phase (u' p a u' p a )^, the Reynolds shear stress 
profile of the fluid phase (u' { v' f ) as well as the correlation between the streamwise and wall-normal 
particle velocity fluctuations {u' v v' p ). All components of the fluid velocity fluctuation covariances 
deviate only marginally from the single phase counterparts when normalized both in outer and 
wall scales. This confirms that the presence of the particles, at the presently adopted parameter 
values, has a negligible influence on the one-point statistics of the fluid velocity field. 

Concerning the dispersed phase, in the near wall region within y + < 36, the particle velocity 
fluctuation intensity is observed to be less than that of the fluid fluctuation in the streamwise 
direction. The reverse is true for the cross-stream components, i.e. in the wall-normal and 
spanwise directions particle velocity components are observed to fluctuate more strongly than 
fluid velocities. An exception are those particles which are in the near proximity of the wall 
(with their centers located at wall normal distance within one particle diameter from the wall) 
where (v' 2 ) ((uip)), is also smaller than (v' f 2 ) ((w' f 2 )). Note that the apparent 'jump' from the 
first to the second data point in figure 6(6) is due to the width of the bins, as can be seen from 
the insert where much finer sampling bins where used and smoother curves are obtained adjacent 
to the wall. The insert also demonstrates the limited amount of particle statistics for locations 
y + > 20. The overall trend of more isotropic particle velocity fluctuations as compared to the 
fluid phase can be explained by the action of particle collisions which tend to re-distribute the 
kinetic energy among the components. Furthermore, as already noted by Kaftori et al. (1995b), 
the action of gravity in the negative wall-normal direction leads to an additional enhancement 
of the fluctuation intensity of the wall-normal particle velocity component. On the other hand, 
the profiles in figure 6(b) (particularly the insert) demonstrate the effect of the finite size of the 
particles: the damping of the particle velocity fluctuations very close to the wall is felt at larger 
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distances than what is experienced by the fluid, and the wall-normal component has a nearly 
vanishing fluctuation intensity at center locations around y + = 5. 

Let us turn to the covariance between streamwise and wall-normal velocity fluctuations. 
Figure 7(a) shows that the fluid and the particle phase exhibit similar values in the near-wall 
region y/h < 0.2. Upon closer inspection (figure lb and insert) it can be seen that the particle 
'Reynolds stress' tends towards zero at the wall faster than the fluid counterpart, and a cross- 
over is observed at y + 7, which again manifests the finite-size of the particles. However, we 
detect a clear local maximum of (u' p v' p ) at a wall-distance of y + ps 10. The insert shows that the 
particle 'Reynolds stress' first decreases significantly for y + > 10, then rises again for y + > 15. 
Although the fluctuation intensity of the two components ((u' p u' p ) and (v' v' )) does not change 
much in the interval 10 < y + < 20, their covariance does exhibit the observed 'kink'. This 
means that the correlation between the two components is somehow disturbed in that region. 
Therefore, we have more carefully investigated the particle paths in the corresponding region. It 
turns out that particle-particle contact is responsible for the observed behavior. Figure 8 shows 
the number of detected collisions as a function of wall-distance. In computing the counter, only 
those collisions were considered which take place in a predominantly vertical direction (i.e. where 
the line connecting the centers of the two colliding particles is closer to the vertical than to any 
horizontal axis). By this choice we eliminate those collisions which occur between particles in 
motion predominantly in a wall-parallel plane. As can be seen, these collisions are concentrated in 
the same interval of wall-normal distances as the reduction of particle 'Reynolds stress' observed 
in figure 7(6). We have analyzed the trajectories of all particles colliding in the mentioned region 
(figure omitted), and found that these collisions mostly correspond to either: (i) particles being 
picked up from the wall by high-speed fluid, then colliding at oblique angles with those particles 
being located on their downstream side; (ii) particles approaching the near-wall region from 
the outer flow, then colliding with near-wall particles. In the former case, an ejection (i.e. in 
the second quadrant of the u' , w p -plane) will be reduced in intensity or converted into a third 
quadrant event. In the latter case, the collision will tend to convert what is most probably a 
'sweep' event (fourth quadrant) into a first quadrant event. In both cases, the consequence will 
be a reduction in amplitude of the (overall negative) correlation between u' p and v' p . 

3.1.4 Mean shear of the fluid phase and mean angular velocity of the particle phase 

Figure 9 shows the profiles of the mean fluid shear, d(uf)/dy, and the mean angular particle 
velocity around the spanwise axis, (w P)Z ). The figure suggests that the average particle rotation 
is driven by the mean shear, as the data points are approximately proportional to the mean 
fluid velocity gradient. We observe that (w p , z ) ~ — Ad(u{)/dy with a proportionality factor of 
A s» 0.33 (visual fit). An analogous result has already been observed in particulate flow through 
a vertical channel (Uhlmann, 2008; Garcia- Villalba et al., 2012) where a proportionality factor 
of A m 0.15 was obtained. Note that the sign of the angular particle velocity corresponds to 
forward "rolling" motion. Incidentally, the angular velocity of particles adjacent to the wall is 
significantly smaller than would be required by a condition of traction on the wall for which 
— (wp. z )i?/(u p ) = 1. In the present simulation, the observed value of — (w PiZ )i?/(u p ) in the 
averaging bin adjacent to the wall is equal to 0.3. Please note that there is no need for particles 
to actually obey a traction condition in the present setup even if they remain in contact with 
the wall, since tangential forces were neglected in our contact model. 

3.1.5 Instantaneous flow field and two-point correlations 

Figure 10 shows the instantaneous value of the streamwise velocity component in a wall-parallel 
plane at a wall-distance of y + = 5. The particulate flow result (figure 10a) can be compared 
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Figure 9: Wall-normal profiles of the mean fluid shear (solid line) and mean value of the spanwisc 
component of particle angular velocity (symbols). The proportionality factor is given as A = 0.33. 

to the result of single-phase flow at otherwise identical parameters provided in figure 106. In 
both cases the flow field exhibits the well-known streaks of high- and low- speed fluid regions. 
From visual observation the result of the two-phase flow differs little from that of the single- 
phase flow and the particles do not seem to alter the turbulence structure much at the given 
wall distance. A more sensitive measure of local flow field modification is provided by the 
dissipation rate of fluctuating kinetic energy, e = 2vs' i ^s' i j, where s'^ — (w-j + u j-.;)/2 * s * ne 
fluctuating rate-of-strain tensor. A map of e is provided in figure 10c for the present case and in 
figure lOd for the corresponding single-phase flow. Note that in the particulate case, where data 
is generated by means of a finite-difference method on a staggered grid, consistent computation 
of the dissipation rate requires interpolation, thereby introducing a certain amount of smoothing. 
Contrarily, the single-phase data is obtained via a pseudo-spectral method in a collocated grid 
arrangement, where dissipation can be directly computed point-wise. Indeed, the figure shows 
that localized contributions to the dissipation rate can be observed in the vicinity of the phase- 
interfaces. However, the largest values of dissipation are in general found in regions outside the 
low-speed streaks, which appear not to be much affected by the presence of particles. Also, note 
that such instantaneous dissipation maps exhibit a considerable temporal variation due to the 
intermittent nature of this quantity. Therefore, a one-on-one comparison is at best qualitative. 
When plane-averaged and averaged in time, the contribution to dissipation stemming from the 
particles should become negligible, since we have seen above that the turbulent kinetic energy is 
practically unaffected by their presence. 

The possible influence of the particles can be further quantified by analyzing the two-point 
auto-correlations of the streamwise velocity fluctuations defined as 

Ru'u'(r x ,y,r z ) = / / / u'(x,y,z)u'(x + r x ,y, z + r z )dxdz (2) 

L x L z {u'u')(y) J Q J Q 

where r x and r z are streamwise and spanwise separations, respectively, and u' refers to a velocity 
fluctuation with respect to the average over wall-parallel planes and time. Note, that in (2) it is 
not distinguished between regions of fluid or solid phase. Figure 11a shows the auto-correlation 
coefficient of the streamwise velocity component as a function of streamwise separations at a wall 
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Figure 10: Visualization of instantaneous flow fields in a wall-parallel plane located at y + — 
5. Graphs (a) and (b) show fluctuations of the streamwise fluid velocity component, (c) and 
(d) show the instantaneous dissipation rate of fluctuating kinetic energy, e = 2vs' i jS / i p where 
= (u'i j +u'j i )/2. Graphs (a), (c) refer to the present particulate flow case (intersections with 
particles are shown in black), and (&), (d) to single phase flow. Note that one matching instant 
in time is shown in (a), (c) and similarly in (6), (d). Moreover, the instant when (a), (c) are 
taken corresponds to the three-dimensional visualization shown in figure 3. 

distance of y + = 10 at zero spanwise separation. In addition the result from the single-phase 
reference case is included in the figure. At increasing separations r x the correlation coefficient in 
the present case is larger than in the single-phase flow and exhibits values of 0.12 at the largest 
streamwise distance compared to a value of 0.04 in the single-phase case. Figure 116 shows the 
correlation coefficient as a function of spanwise separation, r z for zero streamwise separation. 
The magnitude and location of the minima differs, i.e. the minimum is located at r£ = 60 and 
has a value of —0.26 in the present case while in the single-phase flow the minimum is located 
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Figure 11: Graphs (a), (6) show the auto-correlation coefficient of the streamwise velocity com- 
ponent as defined in (2), for separations in the streamwise and spanwise direction, respectively. 
The corresponding one-dimensional energy spectra are shown in graphs (c) and (d) as a function 
of the wavenumber k a (with directional index a = 1 and 3). Solid lines are for the present 
particulate flow case, dashed lines for corresponding single-phase flow. 

at r+ ~ 55 and has a value of —0.19. In contrast to the single-phase flow the auto-correlation 
in the present case reaches positive values for 120 < r+ < 230, leading to a local maximum at 
r+ = 166 (iW= 0.04). 

Since small-scale contributions are not easily observed from the spatial correlation functions, 
we present the corresponding one-dimensional energy spectra in figure ll(c, d). It can be seen 
that particles affect the distribution of kinetic energy among the scales only at scales around 
and below their diameter D. The effect is more pronounced in the streamwise direction, where 
a visible 'bump' is observed. 

In conclusion, a small but distinct effect of the particles on the turbulence structure is found 
in the auto-correlation of the flow field close to the wall. The increase in the magnitude of the 
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Figure 12: The spherical segment surface S over which fluid velocity is averaged in order to 
define a representative instantaneous velocity at the particle location. The spherical segment 
has a radius R s = 3i?. The parallel planes which define the segment are intersecting the sphere 
with radius R s at a distance equal to R on either side of the center, both being parallel to the 
wall-plane. 

correlation coefficient of the streamwise velocity fluctuations for separation in the streamwise as 
well as in the spanwise direction can be explained by a stabilising effect of the particles on the 
near wall structures leading to somewhat more elongated flow structures. 

3.2 Definition of the fluid velocity seen by finite-size particles 

The concept of a fluid velocity "seen" by suspended particles is extensively used in the context 
of particulate flow systems. It represents a simplification in the sense that one supposes particles 
to undergo a certain forcing due to an "incoming" flow, similar to a fixed immersed object in 
a cross-stream, while ignoring the modification to the flow due to the particle itself. When 
considering point particles, particularly in the case of one-way coupling, the characteristic fluid 
velocity is simply taken as the fluid velocity evaluated at the particle's location. 1 Finite-size 
particles however, modify the carrier flow around them by constraining the fluid velocity at the 
fluid/solid interface. Consequently, there is no unique definition of a representative fluid velocity 
seen by finite-size particles. 

Previously, there have been several attempts to define the undisturbed fluid velocity seen by 
finite-size particles. Bagchi and Balachandar (2003) simulated isotropic turbulence swept past a 
fixed sphere and Zeng et al. (2008) considered non-periodic channel flow with a fixed sphere. In 
both of these studies, the undisturbed fluid velocity at a particle location was estimated from the 
fluid velocity at the same location of a separately computed turbulent flow without the sphere. 
Merle ct al. (2005), using Taylor hypothesis, approximated the undisturbed fluid velocity at a 
certain downstream location (2. 51? from the center of the particle) and at an earlier time. Note 
that all of these three studies considered fixed spheres. The authors remark that these approaches 
are only adequate for estimating the velocity seen by particles with sizes smaller than the flow 
scales and they question the applicability in the case of finite-size particles. Lucci et al. (2010) 
considered mobile finite-size particles in decaying isotropic turbulence in the absence of gravity. 
They defined the characteristic fluid velocity as an average over a small spherical cap with center 
located in the direction given by the particle velocity (measured in an inertial reference frame). 

the case of two-way coupling, strictly speaking one would need to consider the velocity at the particle 
location in a companion simulation including all particles except the one under consideration. 
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Figure 13: (a)-(c) Contour plots of instantaneous streamwise velocity fields of an unbounded flow 
past a fixed sphere in the steady axisymmetric and the unsteady wake regimes, (a) i?eoo = 50, 
(6) Reryo = 140, and (c) Re^ — 350. Contourlines correspond to (0:0.1:1) times the magnitude 
of Uf i00 . The circle indicated by a dashed line has a radius of 3R. (d) The computed value uf 
of the axial component of the instantaneous fluid velocity "seen by the particle" , as given by 

the definition (3), plotted as a function of R s . ( ) Reoo = 50; ( ) Re^ = 140; ( ) 

R eoo = 350. 



In turbulent flows, however, this choice of directional bias appears questionable since it is the 
relative (not the absolute) velocity which distinguishes the direction of the "incoming" fluid flow. 

Here, we propose a definition for the characteristic fluid velocity in the vicinity of the particles 
which avoids the directional assumption. The instantaneous fluid velocity uf in the vicinity 
of the ith particle is approximated by the average of the velocity of the fluid located on a 
spherical surface S of radius R s centered at the particle's center location Xp . In order to avoid 
sampling bias due to the inhomogeneity in wall-normal direction in the considered channel flow, 
the averaging spherical surface is trimmed by two wall-parallel planes at one particle radius R 
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below and above the center leading to the definition of a spherical surface segment as shown 
in figure 12. Any point / on this surface has a position vector Xj (i) = Xp (f) + R s n where n 

is an outward pointing unit vector normal to S at xj f\t). With the surface S discretized by 
Ni Lagrangian marker points, the average velocity of the fluid instantaneously located on the 
surface S is defined as (using the indicator function </>f(x, t) defined in 18): 

1 Nl 

uf(xW,i? 8 ,t) = ^ y —^^(xW(t),i)u f (xf ) (t),i) Vi = l...iV p , (3) 
n S (*) i=i 

where the counter rig '(i) is defined as 

JVi 

nf(t)=J2M4 ] '(*)»*)• (4) 
i=i 

Note that the presence of </>f in (3) avoids undesired sampling of velocity data inside a neighboring 
solid particle. Generally x| 1 ' it) does not coincide with the fixed Cartesian grid where the velocity 
field Uf(x, t) is available. Therefore we determine Uf (x|^(£), t) by a trilinear interpolation from 
Uf values at the grid nodes in the fluid domain fif(i). 

The choice of the radius R s needs to meet two requirements. First, it should be chosen 
sufficiently large in order to avoid a possible influence of the particle's own near field upon the 
computed value of the surrounding fluid velocity. Secondly, the value of R s should not be chosen 
too large such that the resulting uf is still of direct relevance to the motion of the ith particle. 

3.2.1 Testing the definition in uniform unbounded flow past a fixed sphere 

In order to define a reasonable value for R s , we have applied definition (3) to resolved numerical 
data of a uniform unbounded flow past a fixed sphere. In this simple flow the relative velocity 
is equal to the incoming fluid velocity Uf j00 . Therefore, this test case allows to estimate the 
smallest value of R s which yields a sufficiently close approximation of the relative velocity, i.e. 
for which uf f=a tif )0o - Here the only parameter is the Reynolds number based on the particle 
diameter and the incoming flow velocity, -Reoo = ui^D /v, for which we have chosen three values, 
namely Reoo = 50, Re^ — 140 (both in the steady axisymmetric wake regime), and Reoo — 350 
in the unsteady vortex shedding regime (see e.g. Johnson and Patel, 1999, for details on the 
structure of wakes behind spheres). Please note that the value of the Reynolds number based on 
particle diameter and the apparent velocity lag (i?ei ag = u\ &% D jv) in the horizontal channel flow 
considered in the main part of this manuscript is in the range of 15 to 20. Returning to the fixed 
sphere in unbounded uniform flow, snapshots of the velocity field in a plane through the center 
of the sphere are depicted in figure 13. It should be mentioned that in this unbounded flow 
configuration trimming the spherical averaging surface is not necessary. However, for reasons of 
consistency with the application of definition (3) to wall-bounded flow, we have performed the 
averaging over the previously defined spherical surface segment (cf. figure 12). 

Figure 13d shows the average value of the instantaneous fluid velocity as defined by (3) while 
varying the value of R s . Note that in the unsteady case (Reoo — 350) an average in time over 
several flow fields was performed. As expected, the value of the fluid velocity defined in (3) 
tends towards the value of the incoming fluid velocity Uf >00 for large values of R s . As i? s tends 
towards R, the surrounding fluid appreciates the presence of the sphere and uf tends to zero due 
to the no-slip condition. Depending on the Reynolds number, the degree to which the computed 
value of the surrounding fluid velocity is affected by the presence of the sphere varies somewhat. 
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Figure 14: (a) Wall-normal profile of the mean streamwise fluid velocity (uf) ( ); the mean 

streamwise particle velocity (u p ) (•); and the mean streamwise fluid velocity in the vicinity of 
the particles (uf) as defined in (3)(o). (6) Wall-normal profile of the average streamwise relative 
velocity between the phases (uf el ) as defined in (5) (A), as well as the apparent velocity lag 
u iag = (v>f) — (u p ) (■). In both figures, outer units are shown in the left and bottom axes and 
wall units in the right and top axes. 



However, it is seen that for all the three cases considered, approximately 90% of the value of the 
incoming flow velocity is obtained at a distance of 3i? from the particle center. Therefore, based 
on the above results, we define the velocity "seen" by the finite-size particles through (3) with 
a choice of R s = 3R. Note that a similar analysis of the fluid velocity field in the vicinity of a 
single fixed particle in forced turbulence has been performed by Naso and Prosperetti (2010). 
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Figure 15: A Voronoi tessellation of a periodic horizontal plane corresponding to the particle 
positions of the snapshot shown in figure 10a. 

3.2.2 The contribution of the mean relative velocity to the apparent velocity lag 

After having found a method to determine the fluid velocity "seen" by a particle as given in (3), 
we can now define an instantaneous relative velocity for the zth particle, viz. 

uf el (i )J R S) t) = uf(x» ) ii B) t)-uW(«) ) (5) 

with R s = 3R. The average over all particles and time of uf and of uf ei can be computed by the 
averaging operator (26) defined in B.3. Note that the values of (uf) and (uf el ) discussed in the 
following stem from 70 instantaneous flow fields. Therefore the statistics are based on a smaller 
number of samples than those computed from data collected during run-time, e.g. the values of 
(itf) and (u p ). Furthermore, we have verified that the following results are not very sensitive to 
a variation of the chosen value of R s . 

Figure 14a shows (uf ) in comparison to (u{) and (u p ) in the wall region y + < 36. The 
figure clearly demonstrates that the average velocity of the fluid in the vicinity of the particles 
(uf ) is significantly smaller than the unconditioned fluid velocity (uf) and that the former is 
comparable in magnitude to the mean particle velocity (u p ) for most of the wall-normal interval. 
Figure 146 which shows (uf ci ) in comparison to u\ ag highlights this finding. At y + rs 5, which 
holds most of the particles, (uf el ) « 0.31u T which is considerably smaller than u lag w 2.10u T . 
At larger wall-normal distances the value of (uf ol ) increases slightly, attaining a maximum value 
of approximately 0.91u r although a clear trend with wall-normal distance is difficult to infer 
due to the decrease of the number of data samples away from the wall (cf. figure 4). However, 
it is beyond doubt that for the considered wall distances (uf el ) is considerably smaller than 
iflag. The relative particle Reynolds number Re vc \ = (uf ci )D/v for the particles closest to the 
wall has a value i?e rc i ~ 2.2 which is much smaller than -Rei ag s» 15. The very small value of 
i?e rc i indicates that there are no substantial particle-induced wakes, which is confirmed by the 
conditionally averaged flow field around the particles discussed in more detail in § 3.4. It can be 
concluded from the present results that the dominant contribution to ui ag stems indeed from the 
preferential location of particles in the low-speed regions of the flow, in agreement with previous 
experimental findings (Kaftori et al., 1995b; Kiger and Pan, 2002). 

3.3 Distribution and motion of particles at the bottom wall 

As a combined effect of gravity and turbulent dynamics, particles in the present study are 
observed to spend most of their time residing near the bottom wall. Occasionally they are 
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Figure 16: Probability density function of (a) the normalized Voronoi cell areas and (b) the ratio 
of the distances between the extreme streamwise and spanwise extents of the Voronoi cells. The 
dashed lines represent corresponding distributions of particle positions generated randomly with 
uniform probability. 

entrained and lifted up by the carrier flow and either return to the bottom or are ejected to 
the outer region of the flow. During their residence time at the wall, the behaviour of particles 
(their streamwise and spanwise motion as well as their spatial distribution) is closely related to 
the dynamics of the near wall coherent flow structures. In this section we study the near-wall 
particle behaviour by analyzing particle positions and velocities. 

It is clearly observable from the instantaneous flow field and from the corresponding particle 
positions shown in figure 10a that particles form elongated structures, resembling streamwise 
aligned chains. Visualization of sequences of such images shows that these particle clusters 



23 



1 



0.8 
0.6 

q 0.4- 

j 

0.2 


-0.2 

20 40 60 80 100 

Figure 17: Lagrangian autocorrelation as a function of time lag r sop for the Voronoi area ( ) 

and the Voronoi cell aspect ratio ( ). 

remain in the form of quasi-streamwise aligned streaky structures which maintain a coherence 
over substantial time scales. In line with previous experimental findings (Kaftori et al., 1995a; 
Nino and Garcia, 1996), these particle streaks are observed to reside preferentially in low-speed 
fluid regions. Animations show that, during their travel downstream, particles remain in these 
relatively quiescent low-speed streaks for relatively long time intervals, exhibiting only slight 
spanwise ondulations. In those animations it can also be seen that incoming sweeps of high 
momentum fluid sometimes act to dissolve the particle accumulations, after which those particles 
relatively quickly migrate back to a neighboring low-speed region. 

As a quantitative measure of the spatial distribution of particles, we have carried out a cluster 
analysis based on Voronoi diagrams. Voronoi diagrams are tessellations which partition a space 
based on a set of given center positions (the 'sites') into regions such that each point inside a given 
region is closer to the region's site than to any other site (Okabe, Boots, and Sugihara, 1992). 
Voronoi diagram analysis is widely utilized in many scientific areas and has been introduced as 
a tool for the analysis of preferential particle concentration in turbulent flows (Monchaux et al., 
2010, 2012). 

In order to investigate the spatial distribution of particles in the immediate vicinity of the wall, 
a horizontal plane of the bi-periodic domain is tessellated based on instantaneous streamwise and 
spanwise positions of all particles whose center is located at a wall-normal distance within one 
particle diameter from the bottom wall. Each of the resulting Voronoi regions has an area Ay 
associated to the ith particle. The example given in figure 15 demonstrates that the tesselation 
is area-filling, taking into account the bi-periodicity. The inverse of a Voronoi cell area is directly 
related to the local particle concentration: small areas correspond to high local particle concen- 
tration and large Voronoi cell areas to voids. We have analyzed the distribution of these areas 
compared to the distribution of the Voronoi areas associated with randomly-positioned particle 
sets. The p.d.f. of the Voronoi areas (normalized by the mean value) of randomly positioned 
points is independent of the particle number density (Ferenc and Ncda, 2007). However, finite- 
size particles have the constraint that they cannot overlap, and at close-packing (the maximum 
possible particle number density), all Voronoi areas are identical and have a Dirac delta distri- 
bution. Here we generate reference data from a sufficiently large number of randomly-positioned 
particle sets (with uniform probability in space) applying the constraint that they do not overlap 
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Figure 18: P.d.f.s of the streamwise particle velocity relative to the mean flow Up (t) — (ui) (y p l * (t)) 
conditioned upon the Voronoi cell area Ay* (t) of all particles whose centers are located at a wall- 
normal distance from the bottom Up it) < D. The color coding represents: conditional p.d.f. 

of those particles such that their associated Voronoi cell area Ay\t) < ^vi ( )'■> conditional 

p.d.f. of those particles such that their associated Voronoi cell area 4>(i) >A y ,( ). The 

dashed line corresponds to the unconditioned p.d.f. of the same quantity in this wall-normal 
interval. Note that the p.d.f.s are not centered about their respective mean. 

and adopting a particle number density equal to the one of the present flow case. 

Figure 16(a) shows the p.d.f. of the normalized Voronoi areas computed from 20314 snapshots 
of the particle 'field'. A significant deviation of the p.d.f. from the corresponding randomly- 
positioned particles indicates that the particle distribution at the bottom wall is far from random 
- as already observed visually, e.g. in figure 15. The DNS data exhibits a much higher probability 
of finding very small and very large Voronoi areas than in a random case with uniform probability. 
The fact that there is a high probability for particles in the present case to be associated with 
very small Voronoi cell areas indicates that, consistent with the visual evidence, particles have a 
clear tendency to accumulate. 

In order to provide a measure of the anisotropy of particle accumulation regions, we have 
further analyzed the shape of the Voronoi cells. More specifically, we have examined the slender- 
ness (i.e. the aspect ratio) of these regions by computing the ratio between the largest extent of 
each Voronoi cell in the streamwise direction (l x y) and in the spanwise direction (4.v)- Figure 
16(6) shows the p.d.f. of the aspect ratio of the Voronoi cells along with the p.d.f. of the same 
quantity for the randomly-positioned particles. Again, there is an appreciable difference between 
the two p.d.f.s. On average, the Voronoi regions associated with the particles from the present 
case have a smaller aspect ratio ((Z x ,v/^z,v) = 0.722) when compared to that of the random case 
((£ x ,vAz,v) = 1-062 - which should tend to unity when increasing the sampling size). Moreover, 
there is a higher probability of finding a small aspect ratio than in the random case and vice 
versa. This indicates that the regions are significantly squeezed (stretched) in the streamwise 
(spanwise) directions which means that particles are more likely to be streamwise aligned than 
aligned along the spanwise axis. Therefore, the Voronoi cell aspect ratio analysis confirms the 
visual observation of alignment of particles into streamwise elongated 'streaks'. 

Furthermore, we have tracked the Voronoi area and Voronoi aspect ratio associated to each 



25 




Figure 19: Lagrangian autocorrelation as a function of time lag r sep for the velocity of all particles 
whose centers are located at a wall-normal wall distance $/p (t) < D. The color code indicates: 

( ) streamwise component, a = 1; ( ) wall-normal component, a — 2; ( ) spanwisc 

component, a — 3. The insert shows a closeup of the same data for small separation times. 



particle in time in order to study the behaviour of these quantities in a Lagrangian sense. The 
Lagrangian autocorrelation function of a Voronoi cell quantity as a function of separation time 
T scp can be defined as 

p , r , (Q'v(t)Qy(t + Ts CP )) , . 

^• QvlTsopj ~ (Q' v (t)Q' v (t)) ' {) 

where the quantity Q\r(t) is either the Voronoi area or the Voronoi aspect ratio at time t and 
Q'y(t) = Qv(t) ~ (Qv)- The averaging in (6) is performed over all particles and over time, where 
segments of time signals are used for which the considered particle resides close to the wall, i.e. 

y$\t') <DVt'e [t,t + T sep ] . 

Figure 17 shows the quantity Rl,q v as a function of r scp for both the Voronoi area and 
Voronoi aspect ratio. The data shown in figure 17 demonstrates that the spatial particle distri- 
bution is indeed highly persistent in time. The Voronoi cell area is still significantly correlated 
(with a coefficient value of approximately 0.2) even after 100 bulk time units. The Voronoi cell 
aspect ratio, on the other hand, changes much faster, with the correlation coefficient practically 
vanishing after approximately 40 bulk time units. The Taylor micro scale (given by the oscu- 
lating parabola at zero separation) measures 5.2Tb (1.3Tb) for the cell area (cell aspect ratio). 
The difference in time scales measured for these two quantities becomes obvious from a closer 
inspection of time sequences of Voronoi diagrams. Therein it can be observed that while parti- 
cles typically remain located for a long time inside one specific accumulation region, the shape of 
their associated Voronoi cell fluctuates significantly without necessarily changing the cell's area. 

The intersection points of the two curves in figure 16(a) (DNS data versus random particle 
distribution) can be used to define an objective criterion for determining whether a particle is 
located inside a cluster or inside a void (Monchaux et a!., 2010). For this purpose, the cell 
area corresponding to the lower (Ay) and upper (Ay) intersection points in that diagram are 
considered as threshold values: all particles with associated Voronoi cell area lower than Ay 
are considered to reside in a cluster, those with a volume larger than Ay are related to voids. 
We have computed for each particle the number of times it crosses these threshold values and 
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we have determined the duration of all temporal intervals during which the Voronoi cell area is 
continuously below (above) the lower (upper) threshold value. This type of study allows us to 
estimate the residence time of particles in a region with an 'extreme' concentration (e.g. a cluster 
or a void). It turns out that the average residence time in a cluster (void) is approximately 7.8Tb 
(6.9Tb) while the average frequency of a particle entering a cluster (void) region is 0.06/Tb 
(0.02/Tb). Concerning the Voronoi cell aspect ratio (which exhibits a single intersection point, 
cf. figure 16&), one can define a threshold below which the Voronoi cell's shape can be considered 
as spanwise elongated. The analysis yields an average temporal interval of 5.5Tb during which 
a particle has an associated spanwise elongated Voronoi cell. 

The data from the Voronoi analysis can be further utilized for the purpose of conditionally 
averaging the Lagrangian particle quantities. To this end we have computed the difference 
between the streamwise component of the instantaneous particle velocity and the mean fluid 
velocity at the same wall-normal distance, i.e. itp (t) — (u{)(yp (£)) for the ith particle. The 
p.d.f. of this velocity difference is shown in figure 18 for those particles whose center is located 
within one diameter of the wall plane. The curve for the unconditioned quantity is clearly 
skewed towards positive values (i.e. large positive fluctuations are more probable than large 
negative ones). Conditioning the same quantity upon the Voronoi cell area has a substantial 
effect, which is, however, mostly limited to a different probability for observing positive values. 
In particular, particles located in clusters (voids) have a significantly lower (higher) probability to 
exhibit streamwise velocities which exceed the average fluid velocity. One can conclude from this 
analysis that particles located in void regions have a higher tendency to be located in high-speed 
fluid regions than particles located in accumulation regions. 

Finally, we consider the auto-correlation function of the particle velocity components, which 
is defined as 

("p ,a(*)"p,q( f + Vp )) 

(»;,aW«;,w) 

in analogy to (6). Note that again only those data points are considered which correspond to 
particles remaining near the wall during the entire separation interval (i.e. with center y[ j (£') < 
D \/t' S [t,t + Ts e p]), consistent with the previous discussion. Figure 19 shows that the wall- 
normal particle velocity component decorrelates fastest, followed by the spanwise component 
and (with a much longer time scale) the streamwise component. The corresponding Taylor 
microscales measure 1.41Tb, 0.15Tb and 0.75Tb for the streamwise, wall-normal and spanwise 
components, respectively. The longer time scale of the streamwise particle velocity is owed to the 
anisotropy of the near- wall turbulence, i.e. the persistence of velocity streaks affecting the particle 
motion. A similar observation has been made in vertical particulate channel flow ( Garcia- Villalba 
et al., 2012). Contrary to that configuration, however, the wall- normal component behaves quite 
differently from the spanwise component in the present case. In particular, the curve for i?L,u p ,2 
exhibits significant negative correlations at small separation times, with the first change of sign 
taking place at r scp = 0.43 Tb. This feature is absent in the corresponding correlation curve in 
the case of the vertical particulate channel, and it must therefore be due to the action of gravity 
in the wall-normal direction. More specifically, gravity is counteracting any particle motion away 
from the wall, leading to wall-normal particle velocities which often exhibit a reversal of sign 
after relatively short times. 



Rh,u P A T scp) = — yjt^ — Va = 1,2,3, (7) 



3.4 Particle-conditioned relative velocity field 

In this section we discuss the fluid velocity field conditioned upon the presence of the particles. 
This analysis provides a quantitative measure of the spatial distribution of the low- and high- 
speed fluid regions with respect to particle positions. First, we define the instantaneous relative 
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Figure 20: The average of the particle-conditioned relative streamwise velocity u x (defined in 
equation 10), normalized by ui ag in the near wall region. The solid black lines correspond to a 
value of unity. 



velocity field with respect to the Zth particle as 

u«(x« i)=u^x« *)-u«(t) (8) 

where 

x«(t) = x -x«(i) (9) 

is an instantaneous coordinate relative to the particle's center position. The average of the 
relative velocity field (with respect to the particle center) over time and over all particles located 
within a given wall-parallel slab s is defined as follows: 

u r (x, !/ W) = (u( , )(x( I ),t)> p , 1 (10) 

where (-) p .t is the particle-centered field averaging operator defined in (28). We have carried out 
the averaging procedure by considering all particles whose centers are located at wall-normal 
distances within one particle diameter from the bottom wall. That is, the considered slab has 
a wall normal thickness of Ah — D and its bottom edge coincides with the bottom wall. Note 
that this averaging procedure is identical to the one employed by Garcia- Villalba et al. (2012) 
in a vertically oriented particulate channel flow. 

If the particle distribution at the bottom were homogeneous, there would be no correlation 
between the spatial distribution of the coherent structures and the particle positions. In this 
case, the average of the particle-centered relative velocity component in the streamwise direc- 
tion M r (x, y("}) is expected to tend to the value ui ag (y^), except in the immediate vicinity of 
the particle. Regardless of the particle distribution, the average relative velocity is zero at the 
fluid/solid interface |x| = D/2, as a consequence of the no-slip condition. If, however, parti- 
cles are located at certain preferred regions with respect to the flow structures, or particle arc 
considerably affecting the flow field such that significant wakes exist behind them, then u r is 
expected to be significantly different from ui ag for appreciable regions of the averaging area, 
except at sufficiently large distances where the averaged flow field is expected to be uncorrelated 
with particle positions. We can rule out the existence of significant particle-induced wakes from 
the analysis in § 3.2.2. Thus, locations where u r < ui ag correspond to regions of surrounding 
low-speed fluid and contrarily, locations where u r > ui ag correspond to high-speed fluid regions 
with respect to particle positions. 

Figure 20 shows the quantity u r /ui ag . It is clearly seen that on average particles are located 
in regions where u r /ui ag is below unity. The low-speed fluid region surrounding the particles 
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Figure 21: The same quantity u r /wi a g as in figure 20. (a) shows a streamwise profile through 
the particle center (z = 0), with the insert presenting a zoom around the particle location. (&) 
shows spanwise profiles through the particle center (at x = 0, solid line) and at x — 6h (dashed 
line). Note that the data in (b) is symmetric with respect to z = 0. 

on average is found to be persistent in the streamwise direction, extending for more than 1000 
wall units both upstream and downstream of the particle position. It reaches the extents of 
the averaging plane (limited by the value of half the streamwise period of the computational 
domain L x /2). Laterally however, the low-speed fluid region is seen to be correlated with the 
particle position for only up to about 30 wall units both in the positive and negative spanwise 
directions, beyond which u r /ui ag rises above unity. The general pattern in figure 20 is one of 
spanwise alternating low- and high-speed regions correlated over the entire streamwise length of 
the domain. Figure 21 shows profiles of the same data on axes through the particle center both in 
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Figure 22: Instantaneous velocity fluctuations of the fluid and the particle phase in a wall- 
parallel plane at y + = 5 indicated by the black (fluid) and red (particles) arrows. The colorplot 
represents the streamwise component of fluid velocity fluctuation, normalized by the friction 
velocity (u'f/u T ). Note that backward pointing arrows do not correspond to upstream motion. 
In fact, upstream particle motion has not been observed in the present case. The graph shows 
only a small part of the entire plane, and the fluid velocity vectors have only been plotted at 
every 15th (10th) grid point in the x (z) direction. 

the streamwise and spanwise directions. Note that here the spanwise profile is averaged over the 
positive and negative directions exploiting the symmetry of the flow. At increasing streamwise 
separation x from the reference particle (figure 21a), u r gradually increases from its smallest value 
(in the vicinity of the particle) towards a finite value of approximately 0.7 ui ag at a distance of 
approximately x — 45 D (i + = 330) in the upstream direction and approximately x — 140 D 
(£ + = 1015) in the downstream direction. For larger distances, u r levels off at this value 
and does not approach unity in the present computational domain. The insert in figure 21(a) 
reveals that the average near-field within several particle diameters exhibits some additional 
structure. On the upstream side, a change of slope occurs at a distance of approximately 1.5D 
from the particle center, while the downstream side features a change of sign and an interval of 
negatively- valued u r up to a distance of roughly ID. It should be noted that the asymmetry of 
the profile of the average relative velocity is not the signature of individual particle wakes, since 
the Reynolds number of the flow around the individual particles is much too low (cf. § 3.2.2) and 
the length scales of the asymmetric features are too large. At this point it can only be stated 
that the observed shape of u r on the ai-axis through the particle center might be related to the 
preferential streamwise location of the particles with respect to the flow features of the near-wall 
layer. As exemplified in one snapshot in figure 22, high-speed particles are predominantly found 
at locations near the front of a high-velocity fluid region, thereby contributing to the asymmetry 
of the streamwise profile of u r . This point, however, merits further investigation in the future. 

In the spanwise direction (figure 216), u r increases faster than in the streamwise direction, 
crossing the value u r = ui ag for the first time at a distance of approximately z = 4. 2D (z + = 30) 
indicating the spanwise extent of the low-speed fluid regions. The particle-conditioned relative 
velocity u r attains a maximum value of 1.2ui ag at a distance of approximately z = 8. 2D (z + = 
60); this characterizes the lateral distance at which the neighbouring high-speed fluid region is 
strongly correlated with the particle position. For larger spanwise distances z, u r /ui ag oscillates 
with a small amplitude around unity reflecting the alternating nature of the low- and high-speed 
fluid regions. 
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It can be concluded that the spatial distribution of the coherent structures exhibits a strong 
correlation with particle positions almost throughout the entire streamwise and spanwise extent 
of the computational domain. A computational domain with dimensions much superior to those 
adopted in the present study would be required if a full decorrelation were desired. 



3.5 Conditional average of the local volumetric particle concentration 

In § 3.3 it has been quantitatively shown through Voronoi analysis that particles are inhomoge- 
neously distributed and that they tend to be streamwise aligned. Furthermore, in § 3.4, it has 
been shown that the most probable position of a near-wall particle is inside a low-speed fluid 
region. Here, we investigate the spatial distribution of the local volumetric particle concentration 
conditioned to particle positions. This type of analysis provides additional information on the 
spatial pattern of particle accumulations/voids. 

The indicator function of the solid phase, ijj, is simply related to the fluid phase indicator 
function <pt defined in (18) by the relation 

tf(x,t) = l-&(x,t). (11) 

Its average over wall-parallel planes (using nodal values of the discrete grid) and over time can 
be defined as 

<*><*> -'-s^k- (12) 

where n(yj) is the discrete counter given in (19) with N t replaced by (the number of available 
temporal records). The quantity defined in (12) is equivalent to the average solid volume fraction 
at a given wall- normal location yj of the discrete grid. Note that (ip) (yj) is alternative to the solid 
volume fraction (<f> a }(v^') defined in (25) based on counting particles whose centers are located 
in wall-normal averaging bins. In the limit of infinitcsimally small particles and an infinitely fine 
grid the difference vanishes. For a given finite grid and particle size, however, the definition stated 
here (i.e. in relation 12) yields smoother results, especially when two-dimensional concentration 
fields are considered, as will be done in the following. 

Using the solid phase indicator function (11), one can define an indicator field which is relative 
to the location of the Zth particle, i.e. ip®(x.V\t), where the definition of the relative coordinate 
given in (9) is used. The average of this quantity over time and over all particles located within 
a given wall-parallel slab s corresponds to a particle-conditioned field of the average solid volume 
fraction relative to particle positions. The average is defined as follows 

^x,yM) = (V«(x (i) ,*)>P,* = 1 - . (13) 

Tip 

where (x) and rip^ are the discrete counter field and the sample counter given in (27) and 
(24) respectively (supposing that the temporal records in both definitions are identical, i.e. 
= N^). For a finite number of samples (n p s ^ > 0) we have by definition that ij> £ [0, 1]. 
The particle-conditioned average solid volume fraction ip can also be interpreted as the prob- 
ability that a point x (with respect to the center of a particle) is located inside the solid phase 
domain O p . Note that for monodispersed particles ip is an even function, i.e. "0(x, j/^) — 
V>(— x, y^)- For a random distribution of particles centered within a given slab s, the average 
local solid volume fraction is equal to (ip)(y a ) and has a uniform distribution throughout the 
averaging volume, except for |x| < D/2 (where tp = 1) and in the immediate vicinity due to 
appreciable finite-size effects. Thus ip/ (ip) is closely related to the pairwise distribution function 



31 



(a) x 



200 400 600 800 1000 




10° 10 1 10 2 

x/D,z/D 

Figure 23: (a) Contours of the conditional local particle concentration normalized by the mean 
solid volume fraction (ip / (ip)) averaged over all particles which are located within the near- wall 
region (y p < D), as defined in (13). Only half of the averaging plane is shown due to symmetry. 
Contour levels correspond to values (1:0.5:4). The graph in (6) shows ip(x, 0, 0, y^)/ (ip) (square 
symbol) and rp(0, 0, z, y^)/{ip) (triangular symbol), both with s = 2. The vertical solid line 
indicates the extent of the reference particle's surface, and the vertical dashed line gives the 
distance |x| — D/2 + 2 Ax from the center of the reference particle which marks the limit of the 
range of action of the collisional force. The chain-dotted lines indicate decay rates proportional 
to x~ ' 5 and x^ 15 . 

which gives a measure of the probability of finding a particle at a certain distance away from a 
given reference particle (see e.g. Sundaram and Collins, 1997; Shotorban and Balachandar, 2006; 
Sardina, Schlatter, Brandt, Picano, and Casciola, 2012). 

Figure 23a shows the spatial map of ip/(ip) for the immediate vicinity of the wall (i.e. for 
all particles whose centers are located in the region within y p < D). The shaded part of the 
map shows the regions where the conditional local particle concentration is higher than that of 
a homogenous case (ip/fy) > 1). This means that there is an increased probability of finding 
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a particle in this region when compared to that of a random particle distribution. Contrarily, 
the white region of the map corresponds to The figure clearly shows the elongated particle- 
conditioned preferred accumulation regions which span the entire box in the streamwise direction 
but which have only a small spanwise extent. A repetitive pattern of the regions is clearly 
observable in the spanwise direction forming high conditional particle concentration ridges and 
low conditional particle concentration troughs. 

More quantitative information can be obtained by considering profiles of tp through the 
particle center both in the streamwise and spanwise directions, respectively (figure 236). Note 
that the reference particle's surface is marked by a vertical solid line in the figure, and the 
distance |x| = D/2 + 2Ax below which the particle-particle repulsion force is active (see § 2.1) 
is marked by a vertical dashed line. For increasing distances |x| > D/2 + 2Ax, both in the 
streamwise and spanwise directions, the quantity tp first rapidly increases, takes a maximum value 
of 13. 8(^) (6.8(tp)) in the streamwise (spanwise) direction at a distance of approximately 1.3D, 
and then decays slowly. The decay depends strongly upon the coordinate direction, reflecting 
the anisotropic nature of the pattern as seen in figure 23 a. The fact that the highest values of tp 
are located within the immediate vicinity of the reference particle is consistent with the results 
of the Voronoi area analysis, where we found that particles tend to locally cluster and result in 
much higher than average local particle concentration values. Moreover, the anisotropic spatial 
distribution of the quantity tp, i.e. the fact that the values of tp are larger in the streamwise 
direction than in the spanwise direction for all values of |x|, is again consistent with the picture 
deduced from the analysis of the aspect ratio of Voronoi cells. Furthermore, in the streamwise 
direction the quantity tp persistently remains above the average value (tp) for extended streamwise 
distances exhibiting two regimes. First, a power law regime (tp cx x~@ with /3 = 0.5) for up to 
a distance of about x = 50D (x + = 370), followed by a plateau regime extending for up to 
x = L x /2. This result shows the large-scale nature of the particle clustering in the streamwise 
direction. Contrarily, the quantity tp decays faster in the spanwise direction (with decay rate 
(3 1.5), first crosses the value tp = (tp) at a lateral distace of approximately z = 5.6D (z + « 40) 
and attains a minimum at a distance of approximately z = 9D (z + « 65). It slowly increases 
and attains a second local maximum which is slightly above the value tp = (tp) at a distance of 
roughly z + « 120. The existence of multiple local maxima is a reflection of the ridge-trough 
pattern visible in figure 23 a. 

The decay rate of the radial distribution function with increasing distance from the reference 
particle has been previously taken as a measure of the small-scale clustering intensity of point- 
particles: the higher the value of the decay rate the stronger is the small-scale clustering (see e.g. 
Sardina et al., 2012; Gualtieri, Picano, and Casciola, 2009). Incidentaly, the decay rate of tp in 
the streamwise direction is identical to that of the radial distribution function of point-particles 
in a channel flow configuration (without gravity) , at a comparable value of the Stokes number 
(Sardina et al., 2012) 

3.6 The role of the streamwise vortices on the distribution of particles 

Visual and quantitative evidence provided so far shows that, consistent with previous experi- 
mental findings, particles are inhomogeneously distributed at the bottom wall. This results on 
the one hand in regions of high particle concentrations where particles are observed to form 
streamwise elongated clusters and on the other hand in void regions of very low particle concen- 
tration. It has been shown that particles tend to avoid the high-speed fluid regions and are on 
average residing in the low-speed fluid regions. As mentioned in the introduction, the behaviour 
of particle segregation at the bottom wall is believed to be linked to the dynamics of the quasi- 
streamwise vortices in the near-wall region. These coherent vortical structures are well known 
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Figure 24: The spatial distribution of clockwise and counterclockwise rotating vortical regions 
with respect to the position of a particle in a cross-sectional plane passing through the particle 
center. The vortical regions are educed from A2 X+ (red-colored contours) and A^" - (blue-colored 
contours) as defined in (15). (a) shows the probability distribution conditioned to all particles 
whose centers are located within one particle diameter from the bottom wall, (b) shows the 
same quantity as in (a) but additionally selecting only particles with positive spanwise velocity 
w p > 0. (c) shows the same quantity as in (a) but additionally selecting only particles with 
negative spanwise velocity w p < 0. Note that the positive x-direction points into the paper. 
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Figure 25: An instantaneous fluid velocity field (indicated by black arrows) around a particle 
in spanwise motion (red arrow) in a cross-sectional plane. The background color indicates the 
streamwise vorticity normalized in wall units, i.e. Wf x - The graph shows only a small part of 
the entire plane, and the fluid velocity vectors have only been plotted at every 5th grid point in 
both coordinate directions. 



to dominate the buffer region in boundary-layer-type flows and are responsible for most of the 
turbulence activity therein (Robinson, 199f ; Jimenez and Pinelli, 1999; Adrian, Meinhart, and 
Tomkins, 2000; Schoppa and Hussain, 2000; Adrian, 2007). 

With the aim to further analyze the role of the streamwise vortices with respect to particle 
segregation, we determine the correlation between particle positions and the position of vortical 
regions in their vicinity. For the identification of vortical regions we have adopted the definition of 
vortical structures proposed by Jcong and Hussain (1995). A vortex is defined as a region where 
the second largest eigenvalue A2 of the tensor S 2 + Vt 2 is negative, S and Jl being the symmetric 
and antisymmetric parts of the fluid velocity gradient tensor Vuf. It is well known that in 
the buffer region of wall-bounded flows, where the quasi-streamwise vortices are dominant, the 
streamwise vorticity field u)f >x — d y W[ — d z V{ is highly correlated with the vortical regions, (Kim 
et al., 1987). Here the sign of CJf, x (positive/negative) is used to determine the sense of rotation 
(clockwise/counterclockwise) of the vortices educed by means of A2. Therefore, we define two 
real- valued fields A2 X+ and A^* - which indicate whether a given point lies inside a vortical region 
with, respectively, clockwise or counterclockwise rotation around the streamwise axis, having the 
amplitude |A2(x, f)|, viz. 



\<g*+(x t) = \ ' I "2 1.x, •>)] u.A2\,x.,i) ^ u, u; f>xl x,t,j ^ u ^ 14 ^ 



f (x,f)|A 2 (x,t)| if A 2 (x,i) < 0, Wf, x (x,i) > 
else 

X"*-( x t) = J <M X >*)I-Mx,t)| if A 2 (x,t) < 0, ,w f , x (x,t) < , 14b , 
2 * ' ' [0 else . 

The presence of the fluid indicator function <f>{ in the above definition guarantees that only points 
in the fluid region are selected. Note that we have additionally excluded a spherical region with 
radius 3i? around each particle center in order to consider only the part of the fluid domain 
Of(t) which is not directly affected by the particle's own near-field (cf. § 3.2). Analogous to 
the notation in previous sections, the fields defined in (14) conditioned upon the Ith particle 
are denoted as Ag* ^(x^, t) and A^* - (xO,t). The average of the quantities A^ 4 ^'' and 
^x-,(0 Qver a jj p ar ticles and time (as in 10) is denoted as 

A^ + (x,^)) = (A^ + ' (i) (x«,i))p,t, (15a) 
Ar(i,yW)HAr' (,) (i (i U)) P ,t- (15b) 
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The fields A2 X+ and A2 ~ correspond to averages of the magnitude of (negative- valued) A2, 
conditioned to particle centers located in the wall-normal slab s, further distinguishing between 
positive and negative streamwise rotation of the fluid. Therefore, they give an indication of 
the intensity of vortical structures (with the respective sense of streamwise rotation) on average 
found at a given position with respect to a particle. 

Figure 24(a) shows the two fields defined in (15) in a cross-sectional plane through the center 
of the reference particle, selecting by the choice of the slab s those particles with y p < D. The two 
graphs clearly show that those particles are on average found at spanwise positions in between 
strong vortical regions of opposite sign. The signatures of these prominent counter-rotating 
vortical regions visible in figure 24(a) are located at wall-normal distances of roughly y + = 20 
to 30 measured from the reference particle center. The sense of streamwise rotation of the two 
intense vortical regions is such that a region of low streamwise velocity is induced at around the 
particle location, in agreement with the results of § 3.4. This result is also consistent with the 
common occurrence of streamwise vortices flanking low-speed velocity streaks with a positive or 
negative spanwise shift according to their sign. 

In order to analyze the correlation between lateral particle motion and the presence of coher- 
ent structures, we have further conditioned the fields defined in (14) with respect to the sign of 
the spanwise particle velocity w p . The result is shown in figure 24(6) and 24(c) for positive and 
negative values of w p , respectively. The first observation from these figures is that the condition 
on the sign of spanwise particle velocity is statistically relevant since the resulting fields are 
clearly distinct from those in figure 24(a). Most prominently, we find regions of large values of 
A2 ~ (^2 X+ ) directly above particles for which w p > (w p < 0). The oppositely signed fields 
(A^ 4 " m case of w p > 0, and A^" - for w p < 0) do not exhibit large values in the direct vicinity 
of the reference particle. These results confirm that the spanwise motion of near- wall particles is 
statistically correlated with coherent vortices found in the buffer layer. The data suggests that 
the lateral particle migration towards the low-speed region is indeed a result of the velocity field 
induced by the quasi-streamwise vortices. 

An example of a particle located below a counter-clockwise rotating vortex and migrating 
in the positive spanwise direction can be found in figure 25. This particle is immersed in the 
boundary layer induced by the streamwise vortex immediately located above the particle (around 
y + s» 20, z + 1=3 0). The direction of the depicted particle's spanwise motion is clearly towards 
the location where an ejection of low-momentum fluid takes place at that instant (z + rs 30 in the 
figure), i.e. a low- velocity region. The scenario which is captured in figure 25 is representative 
of many events involving near-bottom particles and quasi-streamwise vortices which we have 
observed in our data-base. 

4 Conclusion 

In the present study we have investigated the motion of heavy spherical particles in horizontal 
open channel flow by means of DNS. The particles are larger than the viscous length scale of 
the near-wall region, requiring full resolution of the flow in their near-field. The gravitational 
acceleration was chosen such that the overwhelming part of the particles is residing near the 
wall-plane. Since the global solid volume fraction in the considered case is very low (5 • 10 -4 ), 
the maximum of the average particle concentration obtained near the wall is still comparatively 
low. The particle-to-fluid density ratio is moderate, leading to a value of the Stokes number 
based on the viscous time scale which is larger than (but of the same order as) unity. 

We have found that the basic statistics of the present fluid velocity field are essentially the 
same as those of single-phase flow at the same Reynolds number. The two-point correlation 
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of the fluid velocity field, however, is somewhat modified, exhibiting slightly larger correlation 
lengths in the near-wall region where the bulk of the particles are located. 

Concerning the dispersed phase, the mean particle velocity is found to be smaller than the 
mean fluid velocity at all wall distances. In previous investigations (Kaftori et al., 1995b; Kiger 
and Pan, 2002) an explanation of this apparent velocity lag based upon preferential particle 
concentration has already been proposed. However, its confirmation hinges on technical points 
such as the definition of a relative inter-phase velocity and the performance of particle-conditioned 
averaging. Therefore, we have investigated the origin of this apparent velocity lag in more detail. 
For this purpose we have first proposed a method to determine a characteristic fluid velocity 
in the vicinity of a particle (i.e. the fluid velocity "seen" by the particle) for the case of finite 
particle sizes. Our definition is based on an average of the fluid velocity over a spherical surface 
(with radius R s ) centered at the particle center. The validity of this definition was tested in the 
case of a single fixed sphere in uniform flow which also allowed us to calibrate the value of R s , 
henceforth set to three times the particle radius. Subsequently, our definition of fluid velocity 
"seen" by the particles was applied to the present horizontal channel flow data. It was found that 
the two phases are instantaneously much closer to equilibrium than indicated by the apparent 
velocity lag, as already observed by previous authors in similar cases (Kaftori et al., 1995b; Kiger 
and Pan, 2002), albeit using less precise measures of the relative velocity. 

We have analyzed the spatial distribution of near-wall particles by means of Voronoi tes- 
selation in a wall-parallel plane. The analysis shows that particles near the wall are strongly 
accumulating into streamwise elongated structures. From the particle conditioned local volu- 
metric concentration field (closely related to the pairwise distribution function) we deduce that 
these accumulation regions are of a very large scale, with their streamwise extent of the order of 
the current domain size and a spanwise repeating pattern over distances of approximately 120 
wall units. The auto-correlation of Voronoi cell areas was found to decay very slowly (Taylor 
microscale of 5.1 bulk time units), which indicates that the particle accumulation regions are ex- 
tremely stable in time. In fact, their time scale is much larger than the time scale deduced from 
the auto-correlation of particle velocities (the Taylor microscale of the streamwise component 
measures 1.4 bulk time units). 

Furthermore, we have computed the particle-conditioned fluid velocity field (for near-wall 
particles) in a coordinate system relative to the particle centers. The average relative velocity 
is found to exhibit spanwise alternating ridges and troughs (approximately 50 wall units apart) 
which are also essentially spanning the entire streamwise domain size. The principal minimum 
of the average relative velocity is centered upon the reference particle, confirming previous ob- 
servations that the particles are preferentially residing inside the low-speed streaks of the buffer 
layer. 

Finally, we have performed a vortex eduction study (based upon the A2 criterion of Jeong 
and Hussain, 1995) with the purpose of elucidating the role of coherent structures upon the 
spanwise motion of particles located near the wall. The particle-conditioned average of the 
negative part of the A2 field (distinguishing between positive and negative values of the streamwise 
vorticity) clearly shows that near-wall particles are preferentially located at spanwise positions 
in between and beneath counter-rotating vortical regions, their sense of rotation being consistent 
with the quasi-streamwise vortices which are known to flank the typical low-speed streaks in the 
buffer region. When the above particle-conditioned A2 field is additionally conditioned upon the 
spanwise particle velocity, it is observed that the direction of spanwise particle motion is strongly 
correlated with the presence of a strong vortical region found directly above the reference particle 
location. The sign of the streamwise vorticity associated with the vortical region found on average 
above the particles is such that the induced velocity at the particle location is consistent with 
the sign of the spanwise particle velocity. This result suggests that the spanwise particle motion 
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near the wall is indeed caused by the quasi-streamwise vortices in the buffer layer. 

The present study opens several new perspectives. The proposed method to determine a fluid 
velocity "seen" by the particles appears potentially suitable for the investigation of fluid-particle 
interaction in a wide range of flow configurations. It has the advantage of being free from an a 
priori bias towards a specific spatial direction. At the same time the definition can be readily 
evaluated with the aid of flow data generated by means of interface-resolved direct numerical 
simulation, or, possibly, data obtained through high-resolution volumetric measurements at the 
particle scale. 

Concerning the interaction of near-wall turbulence and particles in the present case, the dy- 
namical mechanism leading to the formation (and sustenance) of particle accumulation regions 
in low-speed streaks merits further investigation. In order to provide a "cleaner" and more 
tractable approach to the numerical simulation of such interaction mechanisms it appears ben- 
eficial to consider instead of full turbulence in relatively large domains a smaller "laboratory 
flow" . This can be achieved either by turning to turbulent flow in smaller domain sizes, i.e. the 
minimal flow unit of Jimenez and Moin (1991), or, alternatively, to invariant solutions of the 
Navier-Stokes equations (Kawahara, Uhlmann, and van Veen, 2012). These alternatives will be 
considered in future studies. 

The data presented in this work can be found at the following URL: http://www.ifh.kit. 
edu/dns_data. 
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Appendix A Statistical convergence of the fluid phase data 



Fully developed single phase channel flows are characterized by a linear variation as a function 
of the wall-distance of the plane- and time-averaged total shear stress (r) / pf = ud y (u) — (u'v'}. 
Traditionally this linearity critereon is used to check the convergence of statistics. However, in 
particulate channel flows, particles contribute additional terms to the streamwise momentum 
balance modifying the total shear stress (r). In the context of the immersed boundary method, 
the modified Navier-Stokes equations are solved throughout the entire domain while adding an 
appropriate volume force term f to the equations, which imposes the no-slip condition at the 
location of the particles viz. 



d t u + (u • V)u + Vp = zA7 2 u + f 
V • u =0 



x e n . 



(16) 



where u = (u,v,w) T is the vector of fluid velocities (without distinction between the regions 
occupied by the fluid and the solid), p is the pressure normalized with the fluid density and 
f = (fx, fy, fz) T is the localized volume force term. Integrating the streamwise momentum equa- 
tion over wall-parallel planes and in time (supposing statistical stationarity) , and additionally 
integrating in the wall-normal direction (from the wall plane y' = up to an undetermined 
distance y' — y) yields the following equation: 



(r)/Pf 



{fxW 



1 



(17) 



It can be seen from (17) that there is indeed an additional contribution to the momentum balance 
due to the particle-related forcing term. Thus the modified total shear stress (r)* should vary 
linearly for a fully developed particle laden channel flow. Figure 26 shows the wall-normal profile 
of (t)* for the fluid statistics accumulated excluding the transient period. The near perfect match 
with the linear profile indicates that the flow is in a statistically steady state in the considered 
interval. 
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Figure 26: Wall-normal profile of the plane- and time-averaged total shear stress as defined in 
(17). The dashed line indicates the function (1 — y/h). 
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Appendix B Averaging operations 



B.l Wall-parallel plane and time averaging 

Let us first define an indicator function <j)f(x,t) for the fluid phase which tells us whether a given 
point with a position vector x lies inside f2f(t), the part of the computational domain SI which 
is occupied by fluid at time i, as follows: 

<M X >*) = Q dsc . (18) 



'W^EEE^-H, as) 



Furthermore, a discrete counter of fluid sample points in a wall-parallel plane at a given wall- 
distance yj is defined as: 

Nt n x n z 

m—l i—1 k—1 

where N K and N z are the number of grid nodes in the streamwise and spanwise directions, 
N t is the number of time records in the data set, and a discrete grid position is denoted as 
Xijfe = (xi,yj, Zk) T . Consequently, the ensemble average of an Eulerian quantity £ f of the fluid 
phase over wall-parallel planes and in time, while considering only grid points located in the fluid 
domain, is defined as: 

N t N x N z 
m =l i=l k=l 

In the present study wall-parallel plane averaging of Eulerian quantities was performed during 
runtime from a total of Nt = 26700 flow field snapshots. 

B.2 Instantaneous box averaging 

A global counter of samples in the fluid domain can be defined analogously to (19), viz. 

JVx N y W s 

vf^EEE^^'')' ( 21 ) 

i=l j=l k=l 

where N y is the number of grid nodes in the wall-normal direction. The instantaneous average 
of a quantity £ f over the entire volume fi/ occupied by the fluid is then defined as follows: 

^ N x N y N z 

n xyz {i) i=lj=lk=1 

B.3 Binned averages over particle-related quantities 

Concerning Lagrangian quantities related to particles, the domain was decomposed into discrete 
wall-parallel bins of thickness Ah and averaging was performed over all those particles within 
each bin. Similar to (18), we define an indicator function 4>^ n (y) which tells us whether a given 
wall- normal position y is located inside or outside a particular bin with index j, viz. 

«(.)-{ J L <J-» Ahi » <i&h (23) 
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as well as a sample counter for each bin viz. 



771=1 (=1 



where iV^ is the number of available snapshots of the solid phase. From the sample counter we 
can deduce the average solid volume fraction in each bin, viz. 

„U*) _n3 

<«<»"') = ^>ra- < 25 ' 

Finally, the binned average (over time and the number of particles) of a Lagrangian quantity £ p 
is defined as follows: 

AT<p> N 

(Sp)(y (j) ) - -4 E E^(y P ^nK P (n , (26) 



CO ^ ^ 
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supposing that a finite number of samples has been encountered (rip > 0). Accumulated data 
related to Lagrangian quantities amounts to a total of iV t — 20314 particle fields (more than 
10 7 particle samples). A bin thickness of Ah = D/2 was adopted to evaluate the bin-averages 
unless otherwise stated. The distribution of samples can be deduced from figure 4 which shows 
the wall-normal profile of the average solid volume fraction. 

B.4 Particle-centered averaging 

For performing the particle-conditioned averaging of a field £f(x, t), an averaging volume is 
considered such that the coordinate system within the volume is defined relative to the center 
of each of the considered particles. Applying the indicator functions defined above, we define a 
discrete counter field n( s )(xjjj.) which holds the number of samples obtained through averaging 
at a given discrete location with indices k for a given y-slab with index s, viz. 

W t <r) N p 

771=1 ( = 1 

In equation (27) the symbol t m indicates the time corresponding to the mth snapshot in the 
database (comprising a total of iV t snapshots) and Xyj.(t m ) is a discrete coordinate relative 
to the Ith particle's center position at time t m , as defined in (9). The average of the quantity 
£f over time and over all particles located within the given wall-parallel slab is then defined as 
follows: 



Note that the coordinates x^ fe (£ m ) do not necessarily coincide with the fixed grid used in the 

direct numerical simulation. Thus, values of £ f at X-^ fc (i m ) are determined by a tri-linear inter- 
polation from those values at the neighbouring fixed grid locations. 

For evaluating the average defined by (28) a total amount of = 70 flow field snapshots 
evenly distributed in the considered observation interval were used. 
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